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I  .   INTRODUCTION 

The  emphasis  of  this  thesis  is  two-fold.   The  first  point 
is  to  ga i n  insight  into  the  characteristics  of  a  spectral 
primitive  equation  model.   The  second  is  to  attempt  to  take 
advantage  of  the  spectral  expansions  in  an  initialization 
techn  i  que . 

A  spectral  model  uses  a  Galerkin  formulation  with  the 
spherical  harmonics  as  its  choice  of  basis  functions.   The 
basis  functions  are  e i gen f unct i ons  of  the  Laplacian  operator 
in  spherical  geometry  and,  thus,  are  used  to  advantage  for 
global  integrations.   An  immediate  consequence  of  the  Galerkin 
formulation  is  that  spatial  derivatives  of  a  single  wave  are 
computed  without  the  usual  truncation  error  present  in  finite 
difference  methods.   This  means  that  the  spectral  model 
contains  no  linear  phase  speed  errors. 

The  non-linear  terms  are  computed  via  a  transform  grid 
in  physical  space  (Orszag,  1971).   This  procedure  is  respon- 
sible for  making  the  efficiency  of  the  spectral  model  com- 
parable with  finite  difference  methods.   In  addition,  the 
non-linear  terms  may  be  computed  without  aliasing.   This  pre- 
vents non-linear  instability  and  conserves  invariant  integral 
properties.   Physical  processes  are  also  computed  on  the 
transform  grid. 

The  second  objective  of  this  thesis  is  to  test  a  divergent 
initialization  for  the  spectral  model.   Most  static  and 
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variational  initialization  techniques  of  global  primitive 
equation  models  use  a  constraint  to  balance  the  rotational 
part  of  the  wind  with  the  mass  field  while  setting  the 
divergent  part  to  zero.   Phillips  (I960)  has  shown  that  a 
completely  non-divergent  initialization  cannot  eliminate  the 
gravity  waves  which  would  be  generated  and  his  work  suggests 
that  the  qua s i -geost rop h i c  divergence  could  tend  to  suppress 
these  gravity  wave  modes  for  mid-latitude  synoptic  scales. 
In  addition,  if  the  quas i -geost rop h i c  divergence  had  internal 
modes  close  to  those  which  would  be  developed  by  the  model, 
such  things  as  the  large  scale  precipitation  would  be  more 
realistic  during  the  first  few  hours  of  the  model  forecast. 

It  is  an  aim  of  this  thesis  to  develop  a  diagnostic 
velocity  divergence  to  insert  into  the  fields  of  an  otherwise 
completely  balanced  system  and  to  test  the  subsequent  effect 
on  the  spurious  gravity  wave  noise  and  the  developing  diver- 
gence. 

A  spectral  formulation  of  the  primitive  equations  can 
lend  itself  to  techniques  which  initialize  according  to  scale. 
For  example,  the  q uas i -geost rop h i c  divergence  might  be  com- 
puted only  for  the  synoptic  scales  of  the  model.   Furthermore, 
the  model  conveniently  supplies  diagnostic  information  about 
the  various  sea  I es . 

The  final  goal  of  this  thesis  is  to  develop  a  semi-impli- 
cit time  differencing  scheme  and  test  its  effects  during  the 
first  few  model  hours  after  initialization. 
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I  I  .   MODEL  DESCRIPTION 

The  model  described  here  has  been  developed  by  Rosmond 
(1977)  and  the  details  are  presented  here  for  reference 
purposes.   Similar  formulations  have  been  done  by  Hoskins 
and  Simmons  (1975)  and  by  Bourke  (1974).   The  equations  for 
an  inviscid,  adiabatic  and  hydrostatic  atmosphere  may  be 
wr  i  tten  as 

1.1)  ||  =  -V«(C+f  )V  -  k  •  Vx(RTVq  +  or  |^) 

an  *      -*■         •  av    2    v 

1.2)  |jp  =  k»Vx(£+f)V  -  V-(RTVq  +  0  |I)  -  VZ(<}>  +  ^  ) 
,.3)  |S  =  -?  .  V6  -  a  || 

,.4)  |a  =  _D  _?  .  Vq  -  || 

1.5,  a|i=-RT 

where 

C  -  vorticity  (C  =  k  •  V  x  ^) 

D  -  velocity  divergence  (D  =  V  •  V) 

T  -  temperature 

6  -  potential  temperature 

tt  -  surface  pressure 

V  -  horizontal  velocity  vector 

<J>  -  geopotential  height 

R  -  gas  constant 

Cp-  specific  heat 


12 


f  -  Coriolis  parameter 

a  -  vertical  coordinate    (a  =   p/ir) 

a  -  vertical  velocity        (a  =  da/dt) 

q  -  Jlmr 

The    continuity    equation   (1.4)  may    be    refined    by     integrat- 
ing    in    the    vertical     and     imposing    boundary    conditions    on    G; 

-*■ 
o(0)    -    a(l)    =    0.        Introducing    the    notation    G    =    V»Vq,    we    may 

rewr  i  te    Eq .     (1.4)    as 

1.6)     |§  =  -(5  +  G)    . 
dt 

where  the  over  bar  denotes  a  vertical  average. 

The  vertical  velocity,  a,  may  be  obtained  diagnostical ly 
by  substituting  Eq.  (1.6)  into  Eq.  (1.4)  and  integrating  in 
the  vertical  to  obtain 


,o 
1.7)  a  =    (D  --  G)a  ■  /   (G  +  D)  da  . 


/ 


o 
By  defining  a  horizontally  mean  temperature, 


T  =  T  (a)  +  T?(a,  X,  6,  t) 


and  using  the  following  operator 


a(X,Y)  =   '   |f  +  |Y 
.   2  3X   9u 

i-y 
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We    can    now    expand    the    vector    quantities     in    our    basic 
equations    in    spherical     coordinates    (X,n.)    to    arrive    at 


•8)      If 


1.9) 


3t 


=  -a(A,B) 


a(B,-A)    -V2(r2E  +  <j)  +  RT  q) 


1. 10)     ||       =  -a(Ue,V0)   +  9D  -  a  || 


I  ll)    5sl 


=   -(D  +  G) 


1.12)      ot=  -RT 
3a 


where 


•     SV  RT ' 

=   (C+f)U  +  a  |j  +  ^  cos 

r 


r 


3X 


G     = 


U       3£nir  3£n7r 


E     = 


,      2     3X 
l-y 

2  2 

U     +  V 


3y 


U     = 


2 ( I -y   ) 

=  u  cos  r\/r 


V  =  v  cos  n/r 

X  =    longitude 

r|  =    latitude 

y  =  sin  n 
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Equations  (1.8)  -  (1.12)  are  the  basic  equations  used  in 
the  model.   The  equations  are  represented  spectrally  in  the 
horizontal  and  finite  differenced  in  the  vertical.   The  de- 
pendent variables  are  written  in  terms  of  a  triangularly 
truncated  series  of  spherical  harmonics: 


X  =  Z  X™   Y™ 


where  the  summation  is  a  double  sum  over  m,n  for  J  m  j  <_  M  and 
|m|  <_  I    <_  M .   The  separation  is  such  that  the  coefficients, 
X-,  are  functions  of  time  and  vertical  coordinate  and  the 
spherical  harmonics,  Y. ,  are  horizontal  functions  of  space. 
The  normalization  and  orthogonality  properties  of  the  Y. 
allow  the  coefficients  to  be  obtained  as  follows: 


2tt    +1 


X£  =  <X'  Yl> 


4TT  J  J 


X  Y™    dydX 


The  model  used  for  this  study  has  five  vertical  layers 
staggered  as  in  Fig.  I.   The  even  levels  carry  a  and  the  odd 
levels  carry  all  the  dependent  variables  plus  the  diagnostic 
variables,  U  ,  V  . 

The  finite  difference  scheme  for  the  vertical  advection 
terms  is 


.  dX  I 

a 


9a    2(Aa 


7  [<xk+l    "  V*k+l    +   (xk  "  xk-i>V 


This  formulation  conserves  kinetic  energy. 
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Figure  I.   Staggered  5-!ayer  sigma  coordinate  system 
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Following  Arakawa  and  Lamb  (1976),  the  hydrostatic  equa< 
tion  (1.12)  can  be  finite  differenced  as 


•l3)   *k+i  -V'f'kVi  +  6kV 


where 


,+■    ,*        ,R/C      .  B  ,*    ,-  ,R/C    . 

ak=  '  "(ak/(W   P   3nd  ek  =  (ak+l/(V   P  -  [ 


A  boundary  condition  is  formed  by  integrating  the  hydrostatic 
equation  from  a  =  0  to  I  which  gives 

I 


or 


o  0J 

I     <t>da   -   <J>sfc    =    R     J      T    da 


where   <J)  ,    is  the  -terrain  geopotent  i  a  I  .   This  boundary  con- 
dition is  then  finite  differenced  as 


1.14)     I  4>k(Aak>   -  <f>sfc  =  R  I  Tk(Aak) 
k 


Equations  (1.13)  and  (1.14)  are  then  combined  to  give  the 
matrix  equation 


1.15) 


<J>  =   [CUT  +  <j> 


sfc 


which  is  the  finite  difference  form  of  the  hydrostatic  equa' 
tion  where  the  variables  are  now  column  vectors. 
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The  pressure  tendency  equation  (I.  I  I)  can  be  written  as 


la  =  _ 
at 


Z(Gk  +  Dk) (AaR) 
k 


or  in  matr  i  x  form  as 


1.16) 


|3-  =  -[N]  (G+D) 


Similarly  equation  (1.7)  can  be  written  in  matrix  form  as 


1.17) 


a  =  [Z]  (G+D) 


The  thermodynamic  equation  (I. 10)  is  differenced  consis- 
tent with  the  method  employed  in  the  hydrostatic  equation. 
The  potential  temperature  is  written  in  terms  of  temperature 
and  the  last  term  becomes 


do 


11*,^,  CT,^,  (^ — )R/Cp 


2(Aa.  )   U~k+T  'k+l   £ 
k  a 


RTk  -       - 

+  ^  (G.    -  G  -   D). 

CP  K 


k+l 


v+ vvWr-5   P)] 

V-l 


For  the  purposes  of  the  sem i - i mp I  i c i t  formulation  de- 
scribed later  in  the  paper,  we  separate  the  temperature  into 
its  horizontal  average,  T  ,  and  perturbation,  T'.   The 
matrix,  [Y],  is  defined  such  that 


fcfl  "g 


(ii-,R/cP 


-  T, 


k+ 


2(Aak) 


*       ak 

T.    -(— 

k       * 


R/C 


-) 


\  + 


'k-l 


2(Aak) 


T 


k-l 


a.    =  [Y]a 
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The  vertically  averaged  divergence  can  be  written  5 
[NUD,  so  that 

*  +  * 

-Mo  -  ~  D  =  -[Y]a  -  £1  [N]D  . 
P  P 

Substituting    for    a    from    equation     (1.17)    we    obtain 
*  * 

-[Y]o  -  ^  [N]D  =   -MDGG  -CYXZ:  +51     [N])D 
P  P 

=   -[SDG  +  [Q]D   . 
Finally    the    entire    term    can    be    written    as 


.  I8)     a  ||  =  -[S]G+[Q]D-        '  [*k+|  (T^,  ^^  -  T') 

k  ak+l 

a,   R/C         RT.       RT' 

Vi  p  p 


The  nonlinear  terms  are  computed  using  the  transform 
method  suggested  by  Orszag  (1971).   The  longitudinal  direc- 
tion is  done  using  Fast  Fourier  Transform  and  the  latitudinal 
direction  is  done  with  Gaussian  Quadrature.   The  number  of 
latitudes,  N,  and  longitudes,  M,  satisfy: 

N  >  (3J+I )/2    and    M  >  3J+I  . 

The  number  of  points  are  picked  to  ensure  non-aliased  results 

for  quadratic  terms.   Certain  terms  containing  a  have  triple 

products  and  are  not  computed  alias  free.   Experiments  by 

the  author  and  a  previous  study  (Hoskins  and  Simmons,  1975) 

indicate  that  this  source  of  error  is  negligible  for  the 

length  of  the  integrations  performed  in  this  study. 
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III.   SEMI- IMPLICIT  SCHEME 

This  section  gives  details  of  the  sem i - i mp I  i c i t  time 
differencing  scheme  used  in  some  of  the  experiments  of  this 
study.   Using  the  matrices  defined  in  the  model  description 
section,  equations  (I.9)-(I.I2)  can  be  written 

an    ?       * 

2.1)  |£  +  VZ(<f>'  +  RT  q)  =  FQ 

2.2)  ^   _  [q]D  =  Fj 

2.3)  |H-  +  [Nf|D  =  -G 

2.4)  <j)»  =  [CUT 

where     <^  *  =  dp  —  <|>  , 

Fn  and  FT  represent  the  remaining  terms  in  each  equation 
which  are  not  explicitly  separated  out.   Following  Robert 
et  al.  (1972),  all  the  terms  on  the  left  hand  side  of  equa- 
tions (2. I )— (2.4)  are  time  averaged.   We  define  the  following 
averaging  operator  and  represent  the  local  tendencies  by  a 
Leapfrog  finite  difference  scheme: 

t  =  (   ,t+At   (   ,t-At 


•  «   >  =  «<   ,  -  <   >t+At  -  <   )t-A+ 

at  l  ;  o+^  ;  -     2(A+) 
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Equations  (2. I )— (2.4)  can  then  be  written  as 


2.5)      6+D  =   F     -  V2(?rt  +  RT     q+) 


2.6)      6+T  =   FT  +  [Q]  D+ 


2.7)  6+q   =   -  G  -   [N]   Df 

2.8)  ^"+  =   [Cl  T+ 


Eliminating  ^T   in  equation  (2.5)  by  substitution  of  Eq . 
(2.8),  we  arrive  at  a  time  averaged  set  of  equations: 

2.9)      D+     =   Dt_At  +(At)(FD  -  V2([C]T+  +  RT*  q+)) 

2.10)  Tf     =  T+"At  +(At)(FT  +   [Q]  D+) 

2.11)  q+     =   q+"At  -(At)(G  +  [N]  D1")    . 


Substituting    equations    (2. II)    and     (2.I0)     into    (2.9)    and 
solving    for    D      we    obtain    a    Helmholtz    equation    for    the    time 
averaged    divergence 

2. 1 2)      ([l]+(At)2[B]V2)D+  =   D+_A+  +    (At)[FD-V2([C]T+~A+ 

+(At)[C]FT+RT*q+"A+-(At)    RT*  G] 


where  [l]  is  the  identify  matrix  and  [B]  =  [C][Q]  -  RT  [N], 


For  a  spectral  formulation,  the  equation  is  trivial  to 

solve  since  the  spherical  harmonics  are  e i genf unct i ons  of 

2  -t 

the  Laplacian  V  .   Knowing  D  ,  we  then  compute  the  prognostic 
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variables  T     ,  q      from  equations  (2.6)  -  (2.7).   D 
is  computed  from  the  definition  of  the  time  averaging 
operator : 


Dt+A+  =  2Df  -  D+"4t 


and  the  vorticity,  £     ,     never  entered  the  semi-implicit 
formalism  and  is  computed  from  a  standard  Leapfrog  scheme. 
Calculations  with  the  sem i - i mp I  i c i t  model  required  5%    more 
computer  time  per  time  step.   Stable  integrations  were 
possible  with  a  time  step  of  60  mins. 


22 


IV.   QUAS l-GEOSTROPHIC  DIVERGENCE 

Following  the  analysis  of  Simmons  and  Hoskins  (1976),  a 
diagnostic  divergence   Is  obtained  for  the  case  of  quasi- 
geostrophic  flow  with  the  full  variation  of  the  Coriolis 
parameter  over  the  sphere.   The  relevant  equations  become 


3.1)  S  =  ~  V.  •  V(e+f)  -  *D 


J¥ 


* 


3.2)  *&   -  V2(<j>f  +  RT*q)  =  0 


3'3)  TF=  -%  •  VT  +  Q  (-  a  "^} 

P 


where     V,  -  rotational  part  of  V 

a   -  sigma  dot  at  odd  levels 
/D  =  V  •  (fV(V"2D) ) 


Making  use  of  the  hydrostatic  relationship,  equation 
(1.15),  we  may  write  equation  (3.3)  as 


f-   =  "M(V   •  VT)  +  y   or 


where 


/  .  Cc3  <A   I*  -  •!* 


cp    ° 


To 


2  8 

Taking  V   of  this  equation,  Jf  of  equation  (3.1)  and  -ttt 

of  eq.  (3.2)  and  substituting  we  arrive  at  the  following 

equat  i  on  for  a: 

3.4)  y*a  =  -v"2C  (V.«V(C+f))+  2D>CC](V^-VT)+RtV  ||  . 
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A  second  equation  for  a  is  obtained  from  the  continuity 
equation  in  the  form  of  eq.  (1.17),  which  when  explicitly 
wr  i  tten  out,  i  s 


3.5)  a  =  a(D  +  G) 


/ 


(D  +  G)  da 


where   a  is  specified  on  even  levels. 

The  diagnostic  cr  from  equation  (3.4)  can  be  used  to  com- 
pute the  divergence  using  eq.  (3.5)  provided  we  specify  a 
boundary  condition.   For  the  purposes  of  initialization,  this 
is  most  conveniently  chosen  to  be 


**--<» 


+  G)  E  0 


since  this  will  tend  to  eliminate  external  gravity  waves. 
Applying  this  constraint  to  equations  (3.4)  and  (3.5)  we  ob- 
tain 

3.6)  y*o   =  -V"2(  (V  »V(C+f))+  &D)   +   [CUV.-VT 


7 


3.7)    a  =  - 

o 


The  method  of  solution  will  be  to  compute  cr  using  equa- 
tion (3.6)  from  an  initial  guess  of  the  divergence.   The  even 
level  cr '  s  are  then  found  by  interpolation  and  a  new  guess  of 
the  divergence  is  obtained  by  inverting  equation  (3.7).   The 
updated  divergence  can  be  substituted  back  into  equation  (3.6) 
and  the  procedure  continues  iteratively  until  sufficient 
accuracy  is  obtained. 
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Equations  (3.6)  and  (3.7)  will  now  be  put  into  spectra 
form.   The  equations  contain  three  similar  advective  terms 
It  will  be  sufficient  to  show  the  spectral  form  of  one  of 
them : 

3.8,   V«Pf..^tM»..Ji^.|f^. 


First,  the  terms  are  Individually  transformed  to  physical 
space  by  the  summations: 


J6,rr 


JTl 


,m       I   nimX 


H      V*  ,      ^m  Dm  frImX 


Tk 


l,m 


$  =  ?m 


.      . m  0m  n  i  mX 
lm  ^  P£  I 


,m 


3(C+f) 


m 


8P 
h  8y  "  * 


The  multiplications  in  eq.  (3.8)  are  then  computed  at 

grid  points  and  the  result  is  transformed  back  to  spectral 

space  by  computing  the  inner  product  defined  as: 

+  1  -2tt 


M    I 
3.9)   J,  =  i- 


// 

■1   0 


J(^C+f)Y^  dXdy  E  <J(^C+f),  Y^>  . 


The  integral  is  computed  using  Gaussian  quadrature  in 
the  y  direction  and  Fast  Fourier  Transforms  in  the  X  direc- 
tion.  The  first  term  of  equation  (3.6)  can  then  be 
eva I uated : 

<(V.-V(C+f))  =  V  *(fV(V2J)) 

=  fJ  +  2<Wi  a  (V2j) 

2     8y 
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and  expanded  in  spectral  form  using  the  previous  computation 

,  .  m 
of  J£: 


r  '  +"  r     '  8pm 

<  (VfVCC+f))  Y^>  =  y  y(EJ^)P^y,(l-y2)(ZI^r)  ^i)PM  dy 


where  the  A  integration  has  already  been  done. 
Making  use  of  the  recurrence  relationship 


upM  =  M    M      M   M 


M 


J  4£  - 


2    2 
JT  -  M 


we  can  write  the  first  integrand  as 


MM       MMM       MMM 


and  using  the  orthogonality  of  the  Legendre  polynomials  we 

obtain  for  the  first  integral 
-+l 


/ 


,v    |M  DM,DM    ,  ,M  M        ,M  M 

(Z  J£  P£)PL  dy  =  JL_,    eL+  JL+|   eL+|    . 


The  second  Integral  can  be  done  similarly  using  the 
re  I  at  ionsh  i  p 


3PM 

,t      2,      *  ^^n      M  DM  0    M       DM 

(|_y    )_     =    (£+l,    e£pA_|    .  te£+j    p£+| 


26 


to    obtain 


+  1 


i.l  \) J   Cl-y' 


,M 


3P 


)(Z 


M 

*     DM      A 

P.    dy 


M 


M 


=   J 


-  JM 


l(l+\)    9y        L     M  L+l    L+l  L-l    L 


-I 


Combining    eqs    (3.10)    and    (3.11)    we    obtain 


3.,2>     dftyVlCW.ft  -  ^  e£  4"  ,   ♦  ^  £,   J^ 


The    second    order    operator   X.    D    can    now    be    formulated. 

-MM 
Using    the    notation,    <*D,    Y.  >    =    F       ,    we    can    write    by    analogy 

with    equation    (3.12) 


FM       orL+l      M  nM  L        M        ,M     -, 

FL-  fiC-j-^D^,    +T?reL+|    JL+|]    . 


Thus, 


% ,  m 


and 


■2„     J. 


,M 


L+l      M  ,-M 


L        M        ^M 


<<   D,   YL>  -  <*«D),   YL>  -  ft[— eL  FLH   +ITreL+|    FL+|]   . 


M  M 

Substituting     in    for    F,      .     and    F. . , ,     we    obtain    the    result 
a  L-l  L+l 


3.13)      < 


A.YJJ.  ■ 


n2r.(L+l)L     M     M     >_M 
"  [(I(LTreLeL-l)DL-2 

+    (i^(eM)2+Li+_iL(£M+|)2)DM 


L 


,    ,    L2  +   L 


(L+l) 

M  M 


(L+l (L+2)    "L+l      L+2        L+2 


>  t>"  : . 
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Finally,  eqs  (3.6)  and  (3.7)  may  now  be  fully  expanded 


as 


3.14)  *?  -  ^  +  x&TT^V^"'^  ♦  «V^»3 


u 

3.  15)   ai  =  "  /   (Di  +  GL >  da 


where 


R^  =  [CD  <V^  •  VT,  Y^> 
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V.   DESIGN  OF  THE  EXPERIMENTS 

The  initial  conditions  used  in  the  experiments  of  this 
thesis  were  developed  from  a  ba roc  I  i n i ca I  I  y  unstable  mean 
flow  and  a  small  perturbation  in  zonal  wave  no.  6  (m=6). 
The  mass  field  of  the  mean  flow  was  constrained  to  be  in 
geostrophic  thermal  balance  with  the  mean  wind.   The  vertica 
profile  of  the  mean  wind  was  linear  varying  from  4.5  m/sec 

at  a    =  I  to  49.5  m/sec  at  a    =    0.   The  latitudinal  structure 

2 
varied  as  sin  ( 2ri )  placing  a  jet  max  at  ri  =  45°.   The  verti- 
cal temperature  profile  approximated  the  U.  S.  Standard 
Atmosp  here . 

The  perturbation  in  wave  no.  6  corresponded  to  a  maximum 
v-component  of  .1  m/sec.   The  model  was  integrated  for  192 
hrs  allowing  the  baroclinic  wave  to  grow  to  finite  amplitude 
Fig.  2  shows  a  time  series  of  the  vorticity  during  this 
period.   The  growth  is  strikingly  linear  until  192  hrs  where 
there  is  evidence  of  the  mean  flow  (m  =  0)  changing  due  to 
feedback  from  the  eddies.   The  I  inear  growth  corresponds  to 
a  doubling  time  of  roughly  26  hrs  for  m  =  6.   The  fields  at 
192  hrs  provided  the  balanced  fields  for  the  initialization 
exper  i  men ts . 

Three  types  of  experiments  and  one  control  case  were 
integrated  for  72  hrs  of  model  time. 
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Figure  2.   Vorticity  time  series  during  the  first  192  hrs. 
Vertical  axis  is  logarithmic.   Units  are  sees"' 
Solid  curve  is  the  magnitude  of  wave  no.  0  at 
latitude  60°N.   Dotted  curve  is  m=6,  short 
dashes  are  m=2  and  long  dashes  are  m=l8. 
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The  control  case  consisted  of  initializing  with  the 

balanced  fields  and  integrating.   The  experiments  consisted 

of  the  following: 

Experiment  A:   Initialize  with  the  balanced  fields  but 

no  d  i  vergence 

Experiment  B:   Initialize  with  the  balanced  fields  plus 

the  qua s i -geost rop h i c  divergence 

Experiment  C:   Identical  to  Experiment  A  but  with  a  semi- 
imp  I  i  c  i  t  scheme 

All  experiments  and  the  control  were  first  integrated 

with  a  spectral  truncation  of  wave  no.  21.   These  will  be 

referred  to  as  low  resolution  experiments.   High  resolution 

experiments  were  integrated  for  the  control  and  Experiments 

A  and  B.   These  cases  were  truncated  at  wave  no.  42.   Time 

steps  of  12  mins  for  the  low  resolution  and  6  mins  for  the 

high  were  used.   Experiment  C  was  integrated  with  time  steps 

of  12  and  60  mins.   All  experiments  had  a  weak  Robert  time 

filter  of  0.1  and  no  dissipation  term  was  used. 
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VI  .   RESULTS 

A.   CONTROL  CASE 

Figures  3-6  show  the  initial  fields  at  level  5  for  the 
low  resolution  experiments.   Figures  7-10  are  the  correspond- 
ing high  resolution  fields.   They  are  characterized  by  a 
developing  frontal  zone.   The  advantage  of  the  high  resolu- 
tion is  apparent  in  the  temperature  gradients  across  the 
cold  front  and  warm  front  and  in  the  surface  pressure.   The 
low  in  the  high  resolution  surface  pressure  field  is  5  mb 
deeper.   The  negative  and  positive  areas  of  the  divergence 
and  vorticity  fields  are  approximately  symmetrical  in  the 
low  resolution  fields  and  quite  asymmetric  in  the  high  reso- 
lution fields.   The  region  of  maximum  positive  vorticity 
corresponds  very  well  with  the  maximum  temperature  gradient. 

The  mean  fields  consist  of  an  indirect  cell  induced  by 
the  growing  baroclinic  eddies.  The  indirect  cell  develops 
in  order  to  decrease  the  vertical  wind  shear  which  must  remain 
in  approximate  thermal  balance  with  the  decreasing  pole  to 
equator  temperature  gradient.   The  effect  is  more  pronounced 
in  the  low  resolution  fields.   The  high  resolution  mean 
fields  contain  a  weak  indirect  circulation  and  a  strong 
d  i  rect  cell. 

The  dynamic  situation  during  the  period  of  the  control 
case  is  one  of  a  rapid  cascade  of  energy  to  sma I ler  scales. 
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The  horizontal  scale  of  the  frontal  zone  decreases.   Figures 
11-13  show  this  effect  in  the  high  resolution  temperature 
fields  for  times  out  to  36  hr.   Beyond  this  point,  the 
energy  builds  up  at  the  spectral  limit  of  the  model  since 
there  is  no  dissipation  and  the  72  hr  fields  are  non-meteoro- 
logical (Fig.  14  and  15).   Because  the  energy  blockage  occurs 
at  smaller  scales,  the  high  resolution  case  (Fig.  14)  is, 
paradoxically,  noisier  and  in  greater  error. 

B.   GRAVITY  WAVES 

In  anticipation  of  later  sections,  we  discuss  the 
characteristics  of  inertial  gravity  waves  for  the  model  used 
in  this  study.   The  divergence  is  the  most  sensitive  indica- 
tor of  gravity  wave  motion.   As  discussed  in  Hoskins  and 
Simmons  (1975),  equation  (2.12)  is  the  sem i - i mp I  i c i t  analogue 
of  the  gravity  wave  equation  for  a  multi-layer  model.   The 
eigenvalues  of  the  matrix  [_B2    (eq.  2.12)  give  the  permissible 
e i gen ve I oc i t i es  of  the  gravity  wave  modes.   The  corresponding 
eigenvectors  give  the  vertical  structure.   Figure  16  depicts 
the  eigenmodes  for  the  experiments  in  this  study.   Table  I 
gives  the  periods  of  the  gravity  wave  modes  for  zonal  wave 
no.  6  in  h  rs . 

TABLE  I 


Externa  1 

Fi  rst 

Second 

Thi  rd 

Fou  rth 

3.31 

1  1  .  1 

28.0 

61.7 
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The  total  divergence  may  be  thought  of  as  consisting  of 
a  sum  of  the  individual  eigenmodes: 


4.  I)      D  =  E  ak  Dk 
k 


where  the  D,  are  the  eigenvectors.   The  weighting  coefficients, 
a.,  may  be  obtained  from  the  simultaneous  equations  formed 
by  the  dot  products  of  the  individual  eigenvectors  with  equa- 
tion (4.1).   A  time  series  of  the  coefficients  will  provide 
useful  diagnostic  information  in  later  sections.   Figures  17 
and  18  are  a  time  series  of  the  external  and  first  internal 
modes  for  the  low  resolution  control  case  at  a  point  on  the 
60°  latitude  circle.   It  is  apparent  that  the  coefficients 
contain  both  the  meteorological  modes  and  any  excited  gravity 
modes.   Figure  17  contains  no  high  frequency  oscillation  and 
is  purely  meteorological.   Figure  18  contains  a  meteorologi- 
cal mode  with  a  period  of  roughly  3  days  plus  a  smaller 
amplitude  wave  whose  period  of  approximately  1/2  day  corre- 
sponds well  with  the  calculated  value  for  the  first  internal 
mode  from  Table  I.   The  remaining  internal  gravity  modes 
have  periods  approaching  meteorological  values  and  are  diffi- 
cult to  distinguish. 

C.   EXPERIMENTS  A  AND  B 
I  .   Low  Reso I ut  i  on 

Neglecting  the  divergence  in  the  initialized  fields 
(Experiment  A)  generated  computational  gravity  waves.   Compar- 
ing the  pressure  tendencies  of  the  control  case  (Fig.  19c) 
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TIME 


Figure  I  7. 


Divergence  time  series  of  the  externa 
mode  for  the  low  resolution  control 
case.   Vertical  axis  units  are  sec"1. 
Time  is  in  hrs. 
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Figure  18.   Time  series  of  the  first  internal  mode 
-  otherw  i  se  as  in  Fig.  17. 


50 


48.1 


U.t 


8    jj    g    IS    g    g    s    r   K    Iff    3{ 


M.I 

/L 

24.1 

fl.I 

£ 


iB       £       8       (p       fe 

ft    i   «    Is   s 


IS       f       B       f 

ft     SI     3     Si 


tilt 


. 


46.1 


u 


24. i 


?       IS       if       fe       $       ip       IE 

K     iii     a     3     S     »     S 


?  <F  ? 

w  lu  kJ 

is  $  a 

i*  i*  i* 


i.s 


o 


(0 

o 

O 

en 

- 

~ 

« 

< 

CD 

1 

en 

o 

+- 

ro 

CD 

c 

CJ 

CO 

CD 

E 

— 

iii 

— 

o 

— 

l_ 

l_ 

CD 

+- 

0) 

Q. 

c 

— 

X 

o 

CO 

UJ 

o 

u 
to 

c 

0) 

o 

x: 

— 

XI         •— 

+- 

ro 

ON             +- 

0 

—            3 

T3 

— 

— 

C 

4- 

o 

ro 

l_ 

00 

CD 

CD 

^ 

> 

1_ 

2 

— 

• 

o 

v_* 

co 

— 

u 

CD 

x 

ro 

+- 

c 

L. 

c 

.— 

o 

CD 

•+- 

E 

L0 

— 

— 

+- 

t_ 

<rt> 

CD 

CD 

"v. 

Q-  E 

cr 

X 

•— 

<x> 

LU 

H 

m 

(0         CT* 

o\       — 

CD 

l_ 

U 

O) 

51 


with  those  of  Experiment  A  in  Fig.  19a,  it  is  apparent  that 
high  frequency  oscillations  of  roughly  3  hr  period  are 
excited.   The  initial  magnitudes  of  the  oscillations  are 
approximately  1-2  mb/hr.   The  oscillations  dampen  out  with 
time  as  the  geostrophic  adjustment  takes  place  although  the 
dispersal  of  the  gravity  waves  during  the  adjustment  process 
is  hampered  by  the  longitudinal  periodicity  and  the  symmetry 
of  the  fields  about  the  equator. 

Figure  19b  shows  the  pressure  tendencies  with  the 
qua s i -geostroph i c  divergence  inserted  into  the  initializa- 
tion fields  (Experiment  B).   The  highest  frequency  oscilla- 
tions have  been  reduced  over  Experiment  A  by  a  factor  of 
app  rox  i  mate  I y  2 . 

Comparison  of  the  divergence  modes  in  Figs  20,  21,  22 
allows  us  to  be  more  quantitative.   The  figures  have  the 
control  case  subtracted  out  and  contain  the  gravity  wave 
noise  generated  by  the  approximations  of  each  experiment. 
The  noise  generated  in  the  internal  modes  by  the  non-divergent 
initialization  is  an  order  of  magnitude  larger  than  the  ex- 
ternal mode  and  is  generally  25%    of  the  total  divergence. 
The  figures  show  that  the  external  mode  of  Experiment  B  is 
about  68$  of  Experiment  A.   The  first  internal  mode  shows 
more  noise  in  Experiment  B  initial ly  and  then  decreasing  to 
50$  and  the  second  internal  mode  shows  a  decrease  of  about 
50$  over  Experiment  A.   In  general  the  effect  of  the  quasi- 
geostrophic  divergence  is  to  reduce  the  gravity  wave  noise. 
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2  .   High  Reso I ut  i  on 

Figure  23  shows  the  q uas i -geost roph i c  divergence  at 
level  5  which  was  used  for  Experiment  B.   A  qualitative  com- 
parison with  the  actual  model  divergence  (Fig.  10)  shows  the 
phases  to  be  in  agreement  but  that  the  q uas i -geost rop h i c 
approximation  is  an  overestimate  for  the  perturbation  fields. 
The  zonal  mean  fields  do  not  compare  well.   The  quasi- 
geostrophic  divergence  shows  a  strong  indirect  circulation 
which  is  displaced  poleward  relative  to  the  model  divergence. 
There  is  no  evidence  of  a  direct  cell. 

Figures  24-26  depict  the  noise  generated  in  the 
divergence  modes  for  Experiments  A  and  B.   Generally  they 
show  that  the  gravity  wave  noise  of  Experiment  B  is  66%-12% 
of  that  in  Experiment  A.   The  effect  of  the  quas i -geost roph i c 
divergence  is  to  reduce  the  noise.   However,  this  reduction 
is  less  than  in  the  low  resolution  experiments. 

In  view  of  the  fact  that  the  high  resolution  experi- 
ments were  less  successful  than  the  low  resolution  and  that 
quas i -geostrop h i c  formulations  are  valid  only  for  small 
values  of  the  Rossby  number,  it  is  interesting  to  compare  the 
initial  model  divergence  of  the  control  case  with  the  quasi- 
geostrophic  divergence  in  a  spectral  manner.   Referring  to 
Fig.  27,  it  is  seen  that  there  is  good  agreement  between  the 
model  and  quas i -geost rop h i c  divergences  for  zonal  waves 
corresponding  to  m  =  6,  12,  18.   The  overest i mat i on  by  the 
quas i -geost roph i c  divergence  noted  earlier  is  primarily  in 
wave  no.  6.   The  mean  fields  (m=0)  and  waves  corresponding 
to  m  =  24,30,36,42  do  not  compare  well. 
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One  motive  for  initializing  with  the  quas i -geost ro- 
phic  divergence  is  that,  presumably,  the  "true"  divergence 
would  be  approximated  at  the  start  and  that  the  model  would 
thus  develop  this  divergence  field  earlier  than  a  non- 
divergent  initialization.   Figures  28  and  29  address  this 
point.   They  show  that  at  6  hours  into  the  integrations,  zona 
wave  numbers  6,  12,  18  do,  in  fact,  approach  the  control 
divergence  slightly  faster.   There  is  almost  no  difference 
for  the  higher  wave  numbers  and  the  mean  field  shows  poor 
adjustment  in  either  case.   This  might  be  expected  for  the 
mean  field  since  the  geostrophic  adjustment  process  is  very 
inefficient  at  this  scale  in  a  bounded  region.   On  the  other 
hand,  the  higher  wave  numbers  adjust  very  rapidly  and  after 
6  hrs,  they  are  in  good  agreement  with  the  control  case. 

In  light  of  the  results  showing  good  agreement  in 
the  initial  quas i -geostroph i c  divergence  with  the  control 
divergence  for  zonal  wave  number  6,  12,  18  and  bad  agreement 
for  other  scales,  it  is  natural  to  modify  Experiment  B  such 
that  only  scales  corresponding  to  m  =  6,  12,  18  are  used 
in  the  qua s i -geost rop h i c  divergence.   Figure  30  shows  the 
time  series  of  the  external  divergence  mode  for  this  experi- 
ment.   Figure  31  shows  the  spectral  adjustment  of  the  diver- 
gence at  6  hours.   The  corresponding  results  for  the  unmodi- 
fied Experiment  B  are  Figs.  24  and  29.   The  results  are 
almost  identical  to  the  unmodified  qua s i -geost roph i c  diver- 
gence results. 
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TIME 


Figure  30. 


Time  series  of  the  gravity  wave  noise 
in  the  modified  Experiment  B.   Vertica 
scale  is  sec"1.   Time  is  in  hrs. 
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D.   EXPERIMENT  C 

As  discussed  in  Hoskins  and  Simmons  (1975),  a  linear  sta- 
bility analysis  of  the  sem i - i mp I  i c i t  scheme  shows  that  the 
time  step  is  limited  only  by  the  Rossby  waves.   The  effect 
of  the  scheme  is  to  slow  down  the  faster  moving  gravity  wave 
modes  so  that  the  integrations  do  not  amplify  in  time. 
Slower  moving  atmospheric  motions  are  less  affected.   Table  II 
presents  sample  calculations  of  the  distortion  of  the  periods 
of  the  gravity  waves  present  in  these  experiments  for  two 
time  steps.   These  may  be  compared  with  Table  I.   The  general 
behavior  is  that  the  periods  increase  with  increasing  time 
steps. 

TABLE  I  I 


At 

Externa  1 

First 

Second 

Th  i  rd 

Fou  rt h 

12  m  i  ns 

3.31 

11.14 

28.0 

61.7 

1  54.0 

60  m  i  ns 

5.7 

12.2 

28.5 

61.9 

1  54.  1 

Two  points  will  be  investigated  in  this  experiment.   The 
first  is  that  the  sem i - i mp  I  i c i t  scheme  should  tend  to  smooth 
out  high  frequency  oscillations,  particularly  with  a  large 
time  step.   The  second  is  that  the  sem i - i mp I i c i t  method 
might  possibly  hamper  the  geostrophic  adjustment  process 
since  it  slows  down  the  gravity  waves  which  are  the  mechanism 
for  the  adjustment. 
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Figure  32  shows  the  pressure  tendencies  using  the  semi- 
implicit  method  with  time  steps  of  12  and  60  mins.   Compar- 
ing this  with  the  corresponding  results  for  the  explicit 
integrations  (Fig.  19a),  it  is  evident  that  the  effect  of 
the  sem i - i mp I i c i t  method  is  to  smooth  out  the  3  hr  oscilla- 
tions.  Figure  33  shows  the  relative  effect  of  a  small  vs. 
large  time  step  in  the  sem i - i mp I  i c i t  method  on  the  external 
gravity  waves.   The  larger  period  of  the  60  min  time  step  is 
evident.   The  larger  time  step  initially  amplifies  the  com- 
putational noise,  but  after  24  hrs  it  is  considerably  damped, 
The  smoothing  in  the  pressure  tendencies  is  postulated  to 
come  about  from  the  decreased  velocity  of  the  external 
gravity  mode  assuming  an  advective  time  scale. 

The  second  point  of  this  experiment  is  shown  in  Figs.  34 
and  35.   Comparison  of  the  divergence  coefficients  for  time 
steps  of  12  and  60  mins  shows  no  difference  at  3  hrs.   The 
adjustment  of  the  divergence  to  the  control  case  is  apparent' 
ly  the  same  for  either  time  step. 
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Fi  gure  34. 


3  hr  divergence  coefficients  for  the  low 
resolution  control  case  (dotted  line)  and 
Experiment  C  with  a  [2  m i n  time  step 
(dashed  line).   Vertical  axis  is  latitudina 
wave  no.   Horizontal  scale  is  sec"'. 
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Figure  35.   3  hr  divergence  coefficients  with  a  60 
min  time  step  -  otherwise  as  in  Fig.  34 


72 


V 


SUMMARY  AND  CONCLUSIONS 


The  results  of  this  study  are  based  upon  the  divergence 
developed  by  a  baroclinic  wave  on  a  uniform  earth.   Initiali- 
zation without  the  divergent  wind  component  generated  spurious 
inertial  gravity  waves.   The  internal  state  of  the  initial 
fields  had  nearly  zero  mean  divergence.   Thus,  the  external 
mode   of  the  generated  gravity  waves  was  an  order  of  magni- 
tude less  than  the  total  divergence  and  was  not  significant. 
However,  the  internal  gravity  wave  modes  generated  were  signi- 
ficant.  An  important  source  for  the  external  mode,  namely 
topography,  was  neglected. 

Initializing  with  the  q ua s i -geost rop h i c  divergence 
generally  decreased  the  spurious  gravity  wave  noise  and  aided 
the  ensuing  integrations  to  develop  a  divergence  closer  to 
the  control  case,  but  both  effects  were  sma I  I.   Furthermore, 
the  advantage  of  the  qua s i -geost rop h i c  divergence  decreased 
when  smaller  scales  were  included  in  the  integrations.   In- 
serting the  qua s i -geostroph i c  divergence  into  only  synoptic 
scales  gave  no  advantage.   However,  the  effect  of  this  experi- 
ment may  have  been  reduced  by  the  fact  that  there  were  no 
planetary  scale  waves. 

Using  a  semi-impl  icit  scheme  with  the  same  time  step  as 
an  explicit  integration  had  no  effect  on  the  divergence  modes 
but  did  smooth  out  osci  I lations  in  the  pressure  tendencies. 
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A  larger  time  step  initially  increased  the  amplitude  and 
period  of  the  fast  external  gravity  waves  and  had  less  effect 
on  slower  moving  motions  while  still  stabilizing  pressure 
tendenc  i  es . 

In  general,  the  effect  of  a  divergent  initialization  for 
a  global  primitive  equation  model  is  small.   This  same  con- 
clusion has  been  reached  by  two  similar  studies  -  Houghton 
et  al  (1971)  and  Dey  et  al  (1975).  This  study,  however,  has 
neglected  any  physical  processes  such  as  precipitation  as 
well  as  orography  and  the  effect  of  an  initial  divergence  in 
these  cases  may  well  be  more  important. 
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